function diff = delta_u(m,param)

    % delta(m) in the expansion region

    diff=(1-param.omega_1)./param.alpha./param.c.*m ...
        -(param.omega_0+param.r.*(param.c+param.K)+param.c.*param.zeta)./param.c+param.delta_2.*m.^(1/(1-param.alpha));

    if param.s == 0
        diff=zeros(size(m));
    end

end 